/*
 * GeoManager.cpp
 *
 *  Created on: 14 Oct 2011
 *      Author: allan
 */

#include "GeoManager.h"

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include <Core/Multiphase.h>
#include <Core/Phase.h>
#include <Core/ReactionSystem.h>

#include <Geochemistry/Activity/AqueousActivity.h>
#include <Geochemistry/Activity/GaseousActivity.h>
#include <Geochemistry/Activity/MineralActivity.h>
#include <Geochemistry/GeoDatabase.h>
#include <Geochemistry/GeoEquilibrator.h>
#include <Geochemistry/GeoSystem.h>
#include <Geochemistry/GeoReactor.h>
#include <Geochemistry/Kinetics/MineralKinetics/MineralKinetics.h>

#include <Utils/Algorithms.h>

GeoManager::GeoManager()
{}

void
GeoManager::SetAqueousSpecies(const vector<string>& aqueous_species)
{
	this->aqueous_species = aqueous_species;
}

void
GeoManager::SetGaseousSpecies(const vector<string>& gaseous_species)
{
	this->gaseous_species = gaseous_species;
}

void
GeoManager::SetMineralSpecies(const vector<string>& mineral_species)
{
	this->mineral_species = mineral_species;
}

void
GeoManager::RegisterKinetics(const MineralKineticsInfo& mineral_kinetics_info)
{
	const string mineral = mineral_kinetics_info.GetMineralName();
	
	mineral_kinetics_infos[mineral] = mineral_kinetics_info;
	
	kinetic_species.push_back(mineral);
}

const vector<string>
GeoManager::GetAqueousSpecies() const
{
	return aqueous_species;
}

const vector<string>
GeoManager::GetGaseousSpecies() const
{
	return gaseous_species;
}

const vector<string>
GeoManager::GetMineralSpecies() const
{
	return mineral_species;
}

const Reaction
GeoManager::CreateReaction(const string& reactant, const Multiphase& multiphase) const
{
	auto opt_reaction_data = GeoDatabase::GetInstance().FindReaction(reactant); assert(opt_reaction_data);
	
	// Extract the species and its stoichiometries from the reaction data
	vector<double> stoichiometries;
	vector<string> species;
	
	// Insert the reactant and its stoichiometry (-1.0) into the containers of species and stoichiometries
	stoichiometries.push_back(-1.0);
	species.push_back(opt_reaction_data->reactant);
	
	// Insert the products and their stoichiometry into the containers of species and stoichiometries
	BOOST_FOREACH(auto product, opt_reaction_data->products)
	{
		stoichiometries.push_back(product.first);
		species.push_back(product.second);
	}
	
	// Extract the equilibrium contants and the interpolation temperature points from the reaction data
	vector<double> logK;
	vector<double> T;
	
	BOOST_FOREACH(auto pair, opt_reaction_data->logK)
	{
		logK.push_back(pair.first);
		T.push_back(pair.second);
	}
	
	// Create the reaction equation
	ReactionEquation equation(species, stoichiometries);
	
	// Create the equilibrium constant of the reaction
	EquilibriumConstant K(logK, T);
	
	// Create the reaction
	Reaction reaction(equation, K, multiphase);
	
	return reaction;
}

const Phase
GeoManager::CreateAqueousPhase(const vector<string>& phase_species) const
{
	Phase aqueous_phase("Aqueous", phase_species);
	
	ElectrolyteSolution electrolyte_solution(phase_species);
	
	BOOST_FOREACH(const string& s, phase_species)
		aqueous_phase.SetSpeciesActivity(s, AqueousActivitySetschenowModel(s, electrolyte_solution, 0.1));
	
	BOOST_FOREACH(const string& ion, electrolyte_solution.GetIons())
		aqueous_phase.SetSpeciesActivity(ion, AqueousActivityHKFIonicModel(ion, electrolyte_solution));
	
	aqueous_phase.SetSpeciesActivity("H2O", AqueousActivityHKFWaterModel(electrolyte_solution));
	
	aqueous_phase.SetSpeciesActivity("CO2(aq)", AqueousActivityRumpfCO2Model(electrolyte_solution));
	
	return aqueous_phase;
}

const Phase
GeoManager::CreateGaseousPhase(const vector<string>& phase_species) const
{
	Phase gaseous_phase("Gaseous", phase_species);
	
	BOOST_FOREACH(const string& gas, phase_species)
		gaseous_phase.SetSpeciesActivity(gas, GaseousActivityIdealModel(gas, phase_species));
	
	gaseous_phase.SetSpeciesActivity("CO2(g)", GaseousActivitySpycherReedCO2Model(phase_species));
	
	gaseous_phase.SetSpeciesActivity("H2O(g)", GaseousActivitySpycherReedH2OModel(phase_species));
	
	return gaseous_phase;
}

const Phase
GeoManager::CreateMineralPhase(const vector<string>& phase_species) const
{
	// Build the phase name of the mineral phase (e.g. "Mineral(Calcite|Magnesite|Dolomite|Quartz)")
	string phase_name = "Mineral(";
	
	BOOST_FOREACH(const string& mineral, phase_species)
	{
		phase_name.append(mineral); 
		phase_name.append("|");
	}
	
	*phase_name.rbegin() = ')';
	
	Phase mineral_phase(phase_name, phase_species);
	
	BOOST_FOREACH(const string& mineral, phase_species)
		mineral_phase.SetSpeciesActivity(mineral, MineralActivityIdealModel(mineral, phase_species));
	
	return mineral_phase;
}

const Multiphase
GeoManager::CreateMultiphase() const
{
	vector<Phase> phases;
	
	if(!aqueous_species.empty())
		phases.push_back(CreateAqueousPhase(aqueous_species));
	
	if(!gaseous_species.empty())
		phases.push_back(CreateGaseousPhase(gaseous_species));
	
	BOOST_FOREACH(const string& mineral, mineral_species)
		phases.push_back(CreateMineralPhase({mineral}));
	
	return Multiphase(phases);
}

const ReactionSystem
GeoManager::CreateEquilibriumReactions(const Multiphase& multiphase) const
{
	const vector<string> equilibrium_species = Difference(Filter(multiphase.GetSpecies(), GeoDatabase::IsSecondary), kinetic_species);
	
	vector<Reaction> reactions;
	
	BOOST_FOREACH(const string& reactant, equilibrium_species)
		reactions.push_back(CreateReaction(reactant, multiphase));
	
	return ReactionSystem(reactions);
}

const ReactionSystem
GeoManager::CreateKineticReactions(const Multiphase& multiphase) const
{
	vector<Reaction> reactions;
	
	BOOST_FOREACH(const string& reactant, kinetic_species)
	{
		Reaction reaction = CreateReaction(reactant, multiphase);
		
		if(mineral_kinetics_infos.count(reactant))
		{
			MineralKineticsInfo mkinfo = mineral_kinetics_infos.find(reactant)->second;
			
			MineralKinetics mk(mkinfo, multiphase);
			
			reaction.SetKinetics(mk);
		}
		
		reactions.push_back(reaction);
	}
	
	return ReactionSystem(reactions);
}

const GeoSystem
GeoManager::CreateGeoSystem() const
{
	Multiphase multiphase = CreateMultiphase();
	
	ReactionSystem equilibrium_reactions = CreateEquilibriumReactions(multiphase);
	ReactionSystem kinetic_reactions = CreateKineticReactions(multiphase);
	
	const vector<string> primary_species = 
		Filter(multiphase.GetSpecies(), GeoDatabase::IsBasis);
	
	const vector<string> equilibrium_species = 
		Difference(multiphase.GetSpecies(), Union(primary_species, kinetic_species));
	
	GeoSystem geosystem(multiphase, equilibrium_reactions, kinetic_reactions, 
		primary_species, equilibrium_species, kinetic_species);
	
	return geosystem;
}

const GeoEquilibrator
GeoManager::CreateGeoEquilibrator() const
{
	const GeoSystem system = CreateGeoSystem();
	
	const GeoEquilibrator equilibrator(system.GetMultiphase(), 
		system.GetEquilibriumReactions());
	
	return equilibrator;
}
	
const GeoReactor
GeoManager::CreateGeoReactor() const
{
	const GeoSystem system = CreateGeoSystem();
	
	CannonicalSystem cannonical_system(
		system.GetMultiphase(),
		system.GetEquilibriumReactions(), 
		system.GetKineticReactions(),
		system.GetPrimarySpecies(),
		system.GetEquilibriumSpecies(),
		system.GetKineticSpecies());

	const GeoReactor reactor(cannonical_system);
	
	return reactor;
}

GeoManager::operator GeoSystem() const
{
	return CreateGeoSystem();
}

GeoManager::operator GeoEquilibrator() const
{
	return CreateGeoEquilibrator();
}

GeoManager::operator GeoReactor() const
{
	return CreateGeoReactor();
}
